4.1. Treatment efffect on phenology
#################################
# Summarize annual phenology data
#################################
#get average phenology anomalies per individual
###############################################
#reshape phenology columns to long format
data = melt(pheno, id.vars = c("Species","ID","Treatment","Year"),
variable.name="phenology", value.name="day")
#get means per species and year
phenoMean = data %>%
group_by(Species,phenology,Year) %>%
summarize(Mean = mean(day, na.rm=TRUE))
#merge data
data <-merge(data, phenoMean, all.x=T, by=c("Species","phenology","Year"))
rm(phenoMean)
#get change in phenology relative to mean
data = data %>%
mutate(phenologyAnomaly = ifelse(phenology=="Leaf_out", Mean-day, day-Mean)) %>% #add column
dplyr::select(!c(Mean)) #delete columns
#from long (many rows) to short (multiple columns) format
data3 = dcast(data, Species + ID + Treatment + Year ~ phenology, value.var = c("phenologyAnomaly"))
#summarize years
VariablesVector = c("Leaf_out","Senescence","GSL")
data3 = data.frame(data3 %>%
filter(!Year %in% c("2017")) %>% #delete rows based on condition
group_by(Species,ID,Treatment) %>%
summarize_at(VariablesVector, mean, na.rm = TRUE))
# Add biomass information
#########################
#merge phenology and biomass tables
data3 = merge(data3, dfs$biomass, all.x=T, by=c("Species","ID","Treatment"))
#reshape phenology columns to long format
data3 = melt(data3, id.vars = c("Species","ID","Treatment","Total","Root","Shoot","RSR"),
variable.name="phenology", value.name="phenologyAnomaly")
#get average phenology dates per individual
###########################################
#from long (many rows) to short (multiple columns) format
data1 = dcast(data, Species + ID + Treatment + Year ~ phenology, value.var = c("day"))
#summarize years
VariablesVector = c("Leaf_out","Senescence","GSL")
data1 = data.frame(data1 %>%
filter(!Year %in% c("2017")) %>% #delete rows based on condition
group_by(Species,ID,Treatment) %>%
summarize_at(VariablesVector, mean, na.rm = TRUE))
#merge phenology and biomass tables
data1 = merge(data1, dfs$biomass, all.x=T, by=c("Species","ID","Treatment"))
#reshape phenology columns to long format
data1 = melt(data1, id.vars = c("Species","ID","Treatment","Total","Root","Shoot","RSR"),
variable.name="phenology", value.name="day")
#########################################################################################
#########################################################################################
###########################
# Multivariate linear model
###########################
#Create spring and autumn treatments
####################################
data4 = data3
data4$spring = ifelse(data4$Treatment %in% c("control","autumn"), "no","warming")
data4$autumn = ifelse(data4$Treatment %in% c("control","spring"), "no","warming")
#Multivariate linear model
##########################
resultsLM = data4 %>%
group_by(Species,phenology) %>%
do({model = lm(phenologyAnomaly ~ spring*autumn, data=.) # create your model
data.frame(tidy(model))}) %>% # get coefficient info
filter(!term %in% c("(Intercept)")) %>% # delete rows
dplyr::select(Species, phenology, term, p.value, estimate, std.error) %>% #delete columns
mutate_if(is.numeric, round, 2) %>%
mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) #add significance
#Ordering of factors
term = c("springwarming", "autumnwarming", "springwarming:autumnwarming") #order treatments
resultsLM$term = factor(resultsLM$term, levels=term, ordered=T)
Species = c("Quercus", "Fagus", "Lonicera") #order species
resultsLM$Species = factor(resultsLM$Species, levels=Species, ordered=T)
#order statistics table to match dataframe
resultsLM = resultsLM[with(resultsLM, order(resultsLM$Species, resultsLM$phenology)), ]
data.frame(resultsLM[,-c(7)])
## Species phenology term p.value estimate std.error
## 1 Quercus Leaf_out springwarming 0.00 16.27 1.16
## 2 Quercus Leaf_out autumnwarming 0.06 -2.17 1.09
## 3 Quercus Leaf_out springwarming:autumnwarming 0.79 0.43 1.57
## 4 Quercus Senescence springwarming 0.02 -10.77 4.56
## 5 Quercus Senescence autumnwarming 0.00 13.67 4.31
## 6 Quercus Senescence springwarming:autumnwarming 0.39 -5.43 6.19
## 7 Quercus GSL springwarming 0.27 5.50 4.88
## 8 Quercus GSL autumnwarming 0.02 11.50 4.61
## 9 Quercus GSL springwarming:autumnwarming 0.46 -5.00 6.63
## 10 Fagus Leaf_out springwarming 0.00 13.72 1.66
## 11 Fagus Leaf_out autumnwarming 0.00 -6.73 1.66
## 12 Fagus Leaf_out springwarming:autumnwarming 0.71 0.89 2.37
## 13 Fagus Senescence springwarming 0.12 -3.28 2.05
## 14 Fagus Senescence autumnwarming 0.00 23.11 2.05
## 15 Fagus Senescence springwarming:autumnwarming 0.74 -1.00 2.94
## 16 Fagus GSL springwarming 0.00 10.44 2.21
## 17 Fagus GSL autumnwarming 0.00 16.38 2.21
## 18 Fagus GSL springwarming:autumnwarming 0.97 -0.11 3.16
## 19 Lonicera Leaf_out springwarming 0.00 21.88 1.95
## 20 Lonicera Leaf_out autumnwarming 0.09 -3.54 2.06
## 21 Lonicera Leaf_out springwarming:autumnwarming 0.00 -8.84 2.83
## 22 Lonicera Senescence springwarming 0.04 -4.75 2.19
## 23 Lonicera Senescence autumnwarming 0.00 7.69 2.32
## 24 Lonicera Senescence springwarming:autumnwarming 0.08 5.67 3.19
## 25 Lonicera GSL springwarming 0.00 17.13 2.77
## 26 Lonicera GSL autumnwarming 0.17 4.15 2.93
## 27 Lonicera GSL springwarming:autumnwarming 0.44 -3.17 4.04
#Plot
#####
#create panel labels
dat_text = data.frame(
label = c("a","b","c","d","e","f","g","h","i"),
phenology = rep(c("Leaf_out","Senescence","GSL"),3),
Species = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
term = "springwarming")
#Plot
LMplot = ggplot(data = resultsLM, mapping = aes(x = term, y = estimate,
fill=phenology, alpha=term)) +
geom_bar(position=position_dodge(), stat="identity", color="black") +
geom_errorbar(aes(ymin=estimate-std.error, ymax=estimate+std.error),
width=.2, # Width of the error bars
position=position_dodge(.9))+
geom_hline(yintercept=0) +
geom_text(data = dat_text, mapping = aes(x = Inf, y = Inf,
hjust = 1, vjust = 1,
label = paste("(",label,")",sep=""),
fontface = "bold"))+
stat_summary(geom = 'text', label = resultsLM$sig,
fun.y = mean, vjust = 2) +
scale_alpha_discrete(range = c(1, 0.7)) +
xlab("Phenophase") +
ylab("Effect size (days)") +
scale_fill_manual(values = c("green3", "orange", "blue")) +
facet_grid(Species~phenology) +
plotTheme
LMplot
######################################################
# Check for deviation from control (one-sample T-test)
######################################################
#get means of control group
data2 = data1 %>%
filter(Treatment %in% c("control")) %>%
group_by(Species,phenology) %>%
summarize(Mean = mean(day, na.rm=TRUE))
#merge data
data2 <-merge(data1[!data1$Treatment=="control",], data2, all.x=T, by=c("Species","phenology"))
#get change in phenology relative to control
data2$phenologyChange = data2$day - data2$Mean
#Ordering of factors
Phenology = c("Leaf_out", "Senescence", "GSL") #order treatments
data2$phenology = factor(data2$phenology, levels=Phenology, ordered=T)
Treatments = c("spring", "autumn", "autumn-spring") #order treatments
data2$Treatment = factor(data2$Treatment, levels=Treatments, ordered=T)
Species = c("Quercus", "Fagus", "Lonicera") #order species
data2$Species = factor(data2$Species, levels=Species, ordered=T)
#T-test
resultsTT = data2 %>%
group_by(Species,phenology,Treatment) %>%
do({model = t.test(.$phenologyChange) # create your model
data.frame(tidy(model),
tidy(shapiro.test(x = .$phenologyChange))[1,2])}) %>%
dplyr::select(Species, phenology, Treatment, p.value, estimate, p.value.1) %>% #delete columns
rename(shapiroTest=p.value.1) %>%
mutate(significance = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) %>% #add significance
mutate_if(is.numeric, round, 2)
#order statistics table to match dataframe
resultsTT = resultsTT[with(resultsTT, order(resultsTT$Species, resultsTT$phenology)), ]
data.frame(resultsTT[,-c(7)])
## Species phenology Treatment p.value estimate shapiroTest
## 1 Quercus Leaf_out spring 0.00 -16.27 0.41
## 2 Quercus Leaf_out autumn 0.01 2.17 0.65
## 3 Quercus Leaf_out autumn-spring 0.00 -14.53 0.24
## 4 Quercus Senescence spring 0.00 -10.77 0.36
## 5 Quercus Senescence autumn 0.00 13.67 0.73
## 6 Quercus Senescence autumn-spring 0.52 -2.53 0.37
## 7 Quercus GSL spring 0.03 5.50 0.33
## 8 Quercus GSL autumn 0.01 11.50 0.69
## 9 Quercus GSL autumn-spring 0.01 12.00 0.38
## 10 Fagus Leaf_out spring 0.00 -13.72 0.41
## 11 Fagus Leaf_out autumn 0.00 6.73 0.17
## 12 Fagus Leaf_out autumn-spring 0.00 -7.88 0.77
## 13 Fagus Senescence spring 0.08 -3.28 0.49
## 14 Fagus Senescence autumn 0.00 23.11 0.02
## 15 Fagus Senescence autumn-spring 0.00 18.83 0.22
## 16 Fagus GSL spring 0.00 10.44 0.60
## 17 Fagus GSL autumn 0.00 16.38 0.31
## 18 Fagus GSL autumn-spring 0.00 26.72 0.42
## 19 Lonicera Leaf_out spring 0.00 -21.88 0.08
## 20 Lonicera Leaf_out autumn 0.01 3.54 0.36
## 21 Lonicera Leaf_out autumn-spring 0.00 -9.50 0.98
## 22 Lonicera Senescence spring 0.02 -4.75 0.48
## 23 Lonicera Senescence autumn 0.00 7.69 0.02
## 24 Lonicera Senescence autumn-spring 0.00 8.61 0.01
## 25 Lonicera GSL spring 0.00 17.13 0.44
## 26 Lonicera GSL autumn 0.00 4.15 0.36
## 27 Lonicera GSL autumn-spring 0.00 18.11 0.98
# Plot
######
#create panel labels
dat_text <- data.frame(
label = c("a","b","c","d","e","f","g","h","i"),
phenology = rep(c("Leaf_out","Senescence","GSL"),3),
Species = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
Treatment = "spring")
#Plot
PhenologyChangePlot = ggplot(data = data2, mapping = aes(x = Treatment, y = phenologyChange,
fill=phenology, alpha=Treatment)) +
stat_summary(fun.y = mean,
geom = "bar", color="black") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width=0.2) +
geom_hline(yintercept=0) +
geom_text(data = dat_text, mapping = aes(x = -Inf, y = Inf,
hjust = -0.2, vjust = 1.5,
label = paste("(",label,")",sep=""),
fontface = "bold"))+
stat_summary(geom = 'text', label = resultsTT$significance,
fun.y = mean, vjust = 2) +
scale_alpha_discrete(range = c(1, 0.7)) +
xlab("Treatment") +
ylab("Phenological change (days)") +
scale_fill_manual(values = c("green3", "orange", "blue")) +
facet_grid(Species~phenology) +
plotTheme
PhenologyChangePlot
##########
#Save PDFs
##########
pdf(paste(output.dir,"PhenologyLinearTreatmentModel.pdf",sep="/"), width=7, height=7, useDingbats=FALSE)
LMplot
dev.off()
## quartz_off_screen
## 2
pdf(paste(output.dir,"PhenologyChange.pdf",sep="/"), width=7, height=7, useDingbats=FALSE)
PhenologyChangePlot
dev.off()
## quartz_off_screen
## 2
4.2. Phenology effect on total biomass
##########################
# Univariate linear models
##########################
#reshape biomass columns to long format
data3 = melt(data3, id.vars = c("Species","ID","Treatment","phenology","phenologyAnomaly"),
variable.name="organ", value.name="biomass")
#Transform biomass to percentage anomaly
########################################
#get means per species
data5 = data3 %>%
group_by(Species,organ) %>%
summarize(Mean = mean(biomass, na.rm=TRUE))
#merge data
data3 <-merge(data3, data5, all.x=T, by=c("Species","organ"))
#get precentage change in biomass
data3$relativeBiomass = (data3$biomass / data3$Mean -1) * 100
#order species
Species = c("Quercus", "Fagus", "Lonicera")
data3$Species = factor(data3$Species, levels=Species, ordered=T)
#Extract linear model info
resultsLM = data3 %>%
group_by(Species,phenology,organ) %>%
do({model = lm(relativeBiomass ~ phenologyAnomaly, data=.) # create model
data.frame(tidy(model), # get coefficient info
glance(model),
tidy(shapiro.test(x = residuals(object = model)))[1,2]
)}) %>% # get model info
filter(term != "(Intercept)") %>% # delete intercept info
dplyr::select(Species, phenology, organ, estimate, std.error, p.value, r.squared, p.value.2) %>% # delete columns
rename(shapiro=p.value.2) %>% #rename columns
mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) %>% # add column
mutate_if(is.numeric, round, 2)
#Order phenophases
Organs = c("Total","Shoot","Root","RSR") #order treatments
resultsLM$organ = factor(resultsLM$organ, levels=Organs, ordered=T)
#order statistics table to match dataframe
resultsLM = resultsLM[with(resultsLM, order(resultsLM$Species, resultsLM$organ, resultsLM$phenology)), ]
data.frame(resultsLM)
## Species phenology organ estimate std.error p.value r.squared shapiro sig
## 1 Quercus Leaf_out Total 2.46 0.77 0.00 0.22 0.39 **
## 2 Quercus Senescence Total -2.12 0.49 0.00 0.35 0.95 **
## 3 Quercus GSL Total -1.30 0.67 0.06 0.10 0.35 *
## 4 Quercus Leaf_out Shoot 1.79 0.74 0.02 0.14 0.23 **
## 5 Quercus Senescence Shoot -1.98 0.44 0.00 0.37 0.61 **
## 6 Quercus GSL Shoot -1.53 0.59 0.01 0.16 0.22 **
## 7 Quercus Leaf_out Root 3.03 0.89 0.00 0.25 0.62 **
## 8 Quercus Senescence Root -2.25 0.60 0.00 0.29 0.75 **
## 9 Quercus GSL Root -1.10 0.81 0.18 0.05 0.02
## 10 Quercus Leaf_out RSR 0.98 0.54 0.08 0.09 0.83 *
## 11 Quercus Senescence RSR -0.18 0.39 0.64 0.01 0.16
## 12 Quercus GSL RSR 0.38 0.45 0.40 0.02 0.36
## 13 Fagus Leaf_out Total 0.78 0.51 0.14 0.06 0.28
## 14 Fagus Senescence Total -0.16 0.36 0.66 0.01 0.60
## 15 Fagus GSL Total 0.27 0.41 0.52 0.01 0.44
## 16 Fagus Leaf_out Shoot 0.42 0.52 0.43 0.02 0.31
## 17 Fagus Senescence Shoot 0.09 0.36 0.80 0.00 0.34
## 18 Fagus GSL Shoot 0.37 0.40 0.36 0.02 0.50
## 19 Fagus Leaf_out Root 1.24 0.60 0.05 0.11 0.18 **
## 20 Fagus Senescence Root -0.49 0.43 0.27 0.03 0.57
## 21 Fagus GSL Root 0.13 0.50 0.80 0.00 0.07
## 22 Fagus Leaf_out RSR 0.81 0.42 0.06 0.10 0.28 *
## 23 Fagus Senescence RSR -0.50 0.29 0.10 0.08 0.82
## 24 Fagus GSL RSR -0.14 0.34 0.68 0.00 0.94
## 25 Lonicera Leaf_out Total -0.16 0.20 0.45 0.02 0.39
## 26 Lonicera Senescence Total 0.37 0.30 0.22 0.04 0.44
## 27 Lonicera GSL Total 0.01 0.22 0.96 0.00 0.44
## 28 Lonicera Leaf_out Shoot -0.13 0.22 0.58 0.01 0.75
## 29 Lonicera Senescence Shoot 0.07 0.34 0.84 0.00 0.90
## 30 Lonicera GSL Shoot -0.12 0.25 0.64 0.01 0.79
## 31 Lonicera Leaf_out Root -0.21 0.30 0.50 0.01 0.50
## 32 Lonicera Senescence Root 0.91 0.43 0.04 0.12 0.09 **
## 33 Lonicera GSL Root 0.24 0.33 0.48 0.02 0.18
## 34 Lonicera Leaf_out RSR -0.09 0.33 0.79 0.00 0.27
## 35 Lonicera Senescence RSR 0.87 0.47 0.07 0.09 0.43 *
## 36 Lonicera GSL RSR 0.37 0.36 0.32 0.03 0.26
#Visually inspect model assumptions
data.assumptions = data3
data.assumptions$category = paste(data.assumptions$Species,
data.assumptions$organ,
data.assumptions$phenology, sep="_") #create identifier column
category.list = as.factor(unique(data.assumptions$category)) #create category vector
par(mfrow=c(2,2)) #set plot layout
for (category in category.list){ #loop over categories
tab.subset=data.assumptions[data.assumptions$category==category, ] #create table subset
plot(lm(relativeBiomass ~ phenologyAnomaly, data=tab.subset), main=category)
}
######
#Plots
######
#Effect sizes
#############
#create panel labels
dat_text = data.frame(
label = c("a","b","c","d","e","f","g","h","i","j","k","l"),
organ = rep(c("Total","Shoot","Root","RSR"),3),
Species = c(rep("Quercus",4), rep("Fagus",4), rep("Lonicera",4)),
phenology = "Leaf_out")
#Plot
LMplot = ggplot(data = resultsLM, mapping = aes(x = phenology, y = estimate,
fill=organ, alpha=phenology)) +
geom_bar(position=position_dodge(), stat="identity", color="black") +
geom_errorbar(aes(ymin=estimate-std.error, ymax=estimate+std.error),
width=.2, # Width of the error bars
position=position_dodge(.9))+
geom_hline(yintercept=0) +
geom_text(data = dat_text, mapping = aes(x = -Inf, y = -Inf,
hjust = -0.1, vjust = -1,
label = paste("(",label,")",sep=""),
fontface = "bold"))+
stat_summary(geom = 'text', label = resultsLM$sig,
fun.y = mean, vjust = 2) +
scale_color_viridis(option="viridis",discrete=TRUE) +
scale_alpha_discrete(range = c(1, 0.7)) +
xlab("Phenophase") +
ylab("Effect size (%/day)") +
scale_fill_viridis(option="E",discrete=TRUE) +
facet_grid(Species~organ) +
plotTheme
LMplot
#Biomass
########
#order
Organs = c("Total","Shoot","Root","RSR") #order treatments
data3$organ = factor(data3$organ, levels=Organs, ordered=T)
#create panel labels
dat_text = data.frame(
label = c("a","b","c","d","e","f","g","h","i"),
phenology = rep(c("Leaf_out","Senescence","GSL"),3),
Species = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
organ = "Root")
BiomassPlot = ggplot(data3[data3$organ %in% c("Total","Shoot","Root"),],
aes(x=phenologyAnomaly, y=relativeBiomass, colour=organ)) +
#geom_point(size=0.2) +
geom_hline(yintercept=0) +
geom_vline(xintercept=0) +
geom_smooth(method=lm, fullrange = T, alpha = 0.2) +
labs(x = "Phenology anomaly (days)", y = "Biomass anomaly (%)") +
scale_color_manual(values=c("orange","yellow2","red3"))+
geom_text(data = dat_text, mapping = aes(x = -Inf, y = Inf,
hjust = -0.2, vjust = 1.3,
label = paste("(",label,")",sep=""),
fontface = "bold"), color="black")+
geom_text(data = resultsLM[resultsLM$organ=="Total",],
mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 2.5,
label = paste(round(estimate,1), " %/day", sig, sep="")),
size=3.5, color="orange")+
geom_text(data = resultsLM[resultsLM$organ=="Shoot",],
mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 4.5,
label = paste(round(estimate,1), " %/day", sig, sep="")),
size=3.5, color="yellow2")+
geom_text(data = resultsLM[resultsLM$organ=="Root",],
mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 6.5,
label = paste(round(estimate,1), " %/day", sig, sep="")),
size=3.5, color="red3")+
facet_grid(Species~phenology, scales = "free") +
plotTheme2
BiomassPlot
#RSR
RSRplot = ggplot(data3[data3$organ %in% c("RSR"),], aes(x=phenologyAnomaly, y=biomass)) +
geom_point(size=0.2) +
geom_smooth(method=lm, fullrange = F) +
labs(x = "Day", y = "Root:shoot ratio") +
scale_color_viridis(option="viridis",discrete=TRUE) +
geom_text(data = resultsLM[resultsLM$organ=="RSR",],
mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 2.5,
label = paste("Slope = ", round(estimate,3), " day-1", sig, sep="")),
size=3.5, color="black")+
facet_grid(Species~phenology, scales = "free") +
plotTheme2
RSRplot
#########
#Save PDF
#########
pdf(paste(output.dir,"PhenophasesEffectSizes.pdf",sep="/"), width=8, height=7, useDingbats=FALSE)
LMplot
dev.off()
## quartz_off_screen
## 2
pdf(paste(output.dir,"AveragePhenology.pdf",sep="/"), width=8, height=7, useDingbats=FALSE)
BiomassPlot
RSRplot
dev.off()
## quartz_off_screen
## 2